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Abstract 


We present the first spatially resolved measurements of galaxy properties in the JWST ERO SMACS 0723 field. 
We perform a comprehensive analysis of five 5 «z« 9 galaxies with spectroscopic redshifts from NIRSpec 
observations. We perform spatially resolved spectral energy distribution fitting with BAGPIPES, using six NIRCam 
imaging bands spanning the wavelength range 0.8—5 um. This approach allows us to study the internal structure 
and assembly of the first generations of galaxies. We find clear gradients both in the empirical color maps and in 
most of the estimated physical parameters. We find regions of considerably different specific star formation rates 
across each galaxy, which points to very bursty star formation happening on small scales, not galaxy-wide. The 
integrated light is dominated by these bursty regions, which exhibit strong line emission, with the equivalent width 
of [O irr]2-H reaching up to ~3000-4000 A rest frame. Studying these galaxies in an integrated approach yields 
extremely young inferred ages of the stellar population («10 Myr), which outshine older stellar populations that 
are only distinguishable in the spatially resolved maps. This leads to inferring ~0.5—1 dex lower stellar masses by 
using single-aperture photometry, when compared to resolved analyses. Such systematics would have strong 
implications in the shape and evolution of the stellar mass function at these early times, particularly while samples 
are limited to small numbers of the brightest candidates. Furthermore, the evolved stellar populations revealed in 
this study imply an extended process of early galaxy formation that could otherwise be hidden behind the light of 
the most recently formed stars. 


Unified Astronomy Thesaurus concepts: Extragalactic astronomy (506); High-redshift galaxies (734); Star forming 
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1. Introduction 


By characterizing the physical properties of galaxies in the 
redshift range 5<z< 10, we can study the epoch of 
reionization, when the universe experienced its last phase 
transition (see, e.g., Treu et al. 2013; Mason et al. 2018; 
Robertson 2022, for a review). With well-sampled photometry 
of high-redshift galaxies, we can robustly model their spectral 
energy distributions (SEDs) and infer the properties of their 
stellar populations. Up until now, the rest-frame optical 
emission from galaxies was unavailable at z > 7, having been 
redshifted to the part of the near-infrared spectrum where our 
facilities lacked sensitivity and spatial resolution. While the 
rest-frame UV emission we have had access to is a good tracer 
of unattenuated star formation, it is a poor tracer of stars older 
and less massive than O and B type, which make up the bulk of 
total stellar mass for populations older than a few megayears. 


? Hubble Fellow. 


Original content from this work may be used under the terms 

BY of the Creative Commons Attribution 4.0 licence. Any further 
distribution of this work must maintain attribution to the author(s) and the title 
of the work, journal citation and DOI. 


The latest addition to the space fleet of telescopes, the James 
Webb Space Telescope (JWST), has unprecedented sensitivity 
and spatial resolution in the near-infrared. This has opened up a 
new window into the rest-frame optical emission of high- 
redshift galaxies, allowing us to understand their stellar 
populations for the first time. The Near-Infrared Camera 
(NIRCam; Rieke et al. 2005, 2023) on board JWST allows us 
to reach this spectral range with a unique depth and resolution. 
The Near Infrared Spectrograph (NIRSpec; Jakobsen et al. 
2022) provides high-resolution spectroscopy in the near- 
infrared, which is key to robustly determine the redshift. This 
improves the modeling of the SEDs by constraining a free 
parameter, thus breaking the degeneracies that the redshift has 
with age and dust (see, e.g., Conroy 2013, for a review). 

Lower-redshift studies have been able to resolve galaxies 
and their components (up to z ^ 2 in, e.g., Zibetti et al. 2009; 
Morselli et al. 2019; Nelson et al. 2019; Suess et al. 2019; 
Abdurro'uf & Hirashita 2022; Giménez-Arteaga et al. 2022). 
At higher redshifts, resolved analyses have typically only been 
possible in lensed systems (e.g., Zitrin et al. 2011; Vanzella 
et al. 2017), or in particularly luminous galaxies that break up 
into multiple components (e.g., Matthee et al. 2020; Bowler 
et al. 2022). Nevertheless, integrated photometry has revealed 


THE ASTROPHYSICAL JOURNAL, 948:126 (13pp), 2023 May 10 


the population demographics: stellar mass functions and 
number counts (e.g., Stefanon et al. 2015; Song et al. 2016; 
Stefanon et al. 2021). 

With JWST we can extend for the first time resolved studies 
beyond redshift ~2, introducing the possibility to study in 
unique detail the first generations of galaxies (e.g., Chen et al. 
2023; Hsiao & Coe 2022). These resolved studies will allow us 
to place unique, new constraints on the formation and evolution 
of the first galaxies: their mass assembly histories, modes of 
growth, chemical enrichment, and earliest quenching mechan- 
isms. In order to build a complete picture of galaxy assembly, a 
resolved view of its components is required, to fully understand 
the interplay between the stellar population, dust, and gas in 
z > 6 galaxies. 

The impact of having resolved observations has so far only 
been studied at low redshifts (z C3). Various works have 
compared the inferred physical properties obtained with 
resolved and unresolved observations (see, e.g., Wuyts et al. 
2012; Sorba & Sawicki 2015, 2018; Fetherolf et al. 2020; Vale 
Asari et al. 2020), with diverse conclusions. A resolved 
approach can have multiple advantages, such as decreasing 
degeneracies in the stellar population synthesis models and 
producing more realistic star formation histories (SFHs; Pérez- 
González et al. 2023). In highly star-forming galaxies, the 
outshining of old stellar populations by young ones is of 
particular importance, which a resolved analysis could 
untangle. Sorba & Sawicki (2018) find that the total stellar 
mass derived from integrated SED fitting can be under- 
estimated by a factor of ~5, compared to spatially resolved 
SED modeling, which they attributed to the outshining effect 
by young stars. 

In this paper, we present the observations and first results of 
a multiwavelength analysis of five high-redshift galaxies 
(5 « z « 9) observed with JWST, to study the spatially resolved 
properties of their stellar populations. These galaxies are the 
highest spectroscopically confirmed targets in the JWST ERO 
SMACS 0723 field. Using NIRCam imaging in six bands 
spanning the wavelength range 0.8—5 um, we perform spatially 
resolved SED fitting with BAGPIPES (Carnall et al. 2018). 
There have been multiple works on these targets, albeit always 
from an integrated perspective (e.g., Arellano-Córdova et al. 
2022; Brinchmann 2022; Carnall et al. 2023; Curti et al. 2022; 
Fujimoto et al. 2022; Heintz et al. 2023; Rhoads et al. 2023; 
Schaerer et al. 2022; Tacchella et al. 2022; Trump et al. 2023). 
These papers derive different stellar masses, some of them with 
extremely young ages of the stellar population. Here we present 
the first spatially resolved analysis on these high-redshift 
galaxies, which we propose as a more robust approach to 
accurately calculate their stellar masses. 

This paper is structured as follows. In Section 2, we 
introduce the JWST observations and data reduction procedure. 
Section 3 describes the methodology we use for the modeling 
of the SEDs with BAGPIPES. In Section 4, we present the main 
results and inferred properties of our sample, in both integrated 
and resolved approaches, as well as discussing the implications 
of our analyses. Finally, in Section 5 we present the summary 
and conclusions of our work. Throughout this paper, we 
assume a simplified ACDM cosmology with Hy — 70 km s^! 
Mpc !, 2,,=0.3, and Q4 —0.7. No lensing correction is 
applied throughout this work. Hence, intrinsic stellar masses 
and star formation rates (SFRs) can be obtained by dividing by 
the magnification factor (ji; see Table 1). 
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Table 1 
Redshift, Coordinates, and Magnification Factors for the Five High-redshift 
Galaxies Studied in This Work 


Redshift ID R.A. Decl. Lensing ju 
5.275 8140 110.78804 —73.46179 1.7 
6.383 5144 110.83972 —73.44536 29 
7.663 10612 110.83395 —73.43454 1.6 
7.665 6355 110.84452 —73.43508 2.7 
8.498 4590 110.85933 —73.44916 10.1 


Note. The information is taken from Carnall et al. (2023), where the lensing 
factors (ju) are taken from Oguri (2010). 


2. Data and Observations 


We use the public data of the galaxy cluster SMACS J0723.3 
— 13277] (SMACS 0723), observed by JWST as part of the Early 
Release Observations (ERO; Programme ID 2736; Pontoppidan 
et al. 2022). We use the Near-Infrared Camera (NIRCam; 
Rieke et al. 2005) photometric data from the catalog reduced by 
G. Brammer et al. (2023, in preparation). The data have been 
reduced using the public software package grizli (Brammer 
2019; Brammer & Matharu 2021; Brammer et al. 2022). 
The photometry is corrected for Milky Way extinction 
assuming E(B — V) —0.1909 (Schlafly & Finkbeiner 2011) 
and the Fitzpatrick & Massa (2007) extinction curve. The 
images are point-spread function (PSF) matched to the FA44W 
band on a common 0/04 pixel! scale. We adopt the PSF 
models for use with the grizli mosaics,'? which are based 
on the WebbPSF models. We compute matching kernels for 
each of the PSFs to the F444W PSF using a Richardson —Lucy 
deconvolution algorithm (Richardson 1972; Lucy 1974) and 
then convolve the images with the resulting kernels to match 
the PSF resolution in FA44W. 

In this work we focus on the five galaxies at 5 « z « 9 that 
have spectroscopic redshifts confirmed from NIRSpec observa- 
tions, presented in Carnall et al. (2023). These are detected in 
the six deep NIRCam imaging filters FO90W, F150W, F200W, 
F277W, F356W, and F444W. The targets also have shallower 
imaging data obtained with NIRISS with the F115W and 
F200W filters, which we exclude in the analysis presented in 
this work owing to their lower spatial resolution. The sample 
spans a redshift range from 5.275 for the closest galaxy up to 
8.498 for the highest-redshift object (Carnall et al. 2023). The 
sources have magnification factors between 1.6 and 10.1 in the 
GLAFIC lens models from Oguri (2010). They are not 
significantly distorted by the gravitational lensing so that it 
affects the spatially resolved SED fitting. Table 1 provides the 
basic information for the targets. We do not apply any lensing 
correction to our results, although we report in Table 1 the 
lensing factors for these sources presented in Carnall et al. 
(2023). Figure 1 displays the cutout images of the five galaxies 
in all available observed bands, as well as the color images 
built combining the F150W, F277W, and F444W filters. 


3. Methodology 
3.1. SED Fitting with BAGPIPES 
To model the SED of the individual pixels and derive the 
physical properties, we use the SED fitting code BAGPIPES 


0 https: / [ github.com/gbrammer/grizli-psf-library /tree/main /smacs0723 


THE ASTROPHYSICAL JOURNAL, 948:126 (13pp), 2023 May 10 


F090W F150W F200W 
in 
N 
N 
in i 
II i 
N * 5 
"mE Won 
NS : ae 
II E US aum. 
N á 
O =e n, 
Ne) " 5 : ;' 
N TC j 
II L1 | 
N T ; 
in a " - 
O | i 
‘© P s , 
Es T j 
II EE -e 
N : 
" | a 
ON 5 z 
Y L 
[so] Bo" 
II i 
N 


F277W 


Giménez-Arteaga et al. 


F356W 


GB 
1 
ID8140 
1 kpc 
ID5144 
1 kpc 
ID6355 
1 kpc 
ID4590 


R 
p 


Figure 1. Cutout images of the five SMACS 0723 galaxies in all available NIRCam bands. The cutouts are 1/2 across and centered in the coordinates provided in 
Table 1. The RGB images are built combining the F150W (B), F277W (G), and F444W (R) PSF-matched images (to the F444W band). The physical scale is 


calculated in the lens plane. 


(Carnall et al. 2018). We set the NIRSpec spectroscopic 
redshifts indicated in Table 1, in order to break the degeneracy 
with age and dust that a photometric or uncertain redshift 
would introduce. We use the SPS models by Bruzual & Charlot 
(2003) and include the nebular emission with CLOUDY (Ferland 
et al. 2017), extending the BAGPIPES default grid (which 
normally reaches up to log,y(U) = —2) so that the ionization 
parameter, U, varies in the range —3 < log,,(U) < — 1, since 
z > 6 galaxies display higher ionization parameters than low- 
redshift galaxies (e.g., Curti et al. 2022; Sugahara et al. 2022). 
We assume a Kroupa (2001) initial mass function (IMF) and a 
Calzetti et al. (2000) attenuation curve, in order to reduce the 
number of free parameters in our fits (since other parameter- 
izations such as the Salim et al. 2018 attenuation curve 
introduce extra parameters such as the bump strength and slope 
6). We choose a constant SFH model, following Carnall et al. 
(2023), which seems adequate to fit our galaxies (we obtain 
reduced X? values within the range 0.1-7.5, shown in further 
detail in the following sections). The formation of very young, 
low-mass, and low-metallicity galaxies is likely bursty, and a 
constant SFH accurately resembles this on short timescales. We 
let the maximum age grid vary from 1 Myr to 1 Gyr, to allow 
for the presence of more evolved stellar populations, and 
further limited by the age at the given redshift. We set the 


visual extinction to vary from Ay=0 to Ay=2 and the 
metallicity from 0 to Zo, with uniform priors. Even though the 
metallicity has been calculated with integrated NIRSpec spectra 
for these targets (see, e.g., Brinchmann 2022; Curti et al. 2022; 
Schaerer et al. 2022), they obtain varying results within our 
allowed range. Moreover, we want to allow for spatial variation 
across the galaxy; thus, we do not set Z as a fixed parameter in 
the fit. Finally, we set the lifetime of birth clouds to 10 Myr. 


3.2. Pixel-based Modeling 


In this work we perform SED fitting on a pixel-by-pixel 
basis. With the setup described in the previous subsection, we 
fit the SED and infer the physical parameters of each individual 
pixel. This allows us to recover the 2D distribution of 
properties such as the stellar mass and SFR. We use SExtractor 
(Bertin & Arnouts 1996) to derive the segmentation map of 
each source. The pixels in our maps correspond to physical 
sizes between 180 and 240 pc. In order to fit the SED of 
individual pixels, we impose a signal-to-noise ratio (S/N) 
threshold of 2 on both the F150W and F200W bands, which are 
the noisiest. We find that this threshold is enough to produce 
trustworthy fits, obtaining good reduced y values in the fits of 
individual pixels, as we show in more detail in the following 
section. Pixels that do not fulfill this S/N criterion are not fitted 
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Figure 2. Maps of the galaxy ID8140 at z — 5.275. First column: maps of the empirical colors in AB mag, inferred as the difference of the bands F150W and F277W 
(top) and the bands F356W and F444W (bottom). The contours correspond to the non-PSF-matched F200W image. The physical scale is calculated in the lens plane. 
The FWHM of the F444W PSF is indicated. Rest of the panels: maps of the physical properties inferred with BAGPIPES. The top row, from left to right, shows the 
resulting maps for the SFRD, the logarithm of the SMD, and the visual extinction Ay. The bottom row shows the maps for the mass-weighted stellar age, the UV slope 
(8), and the EW of the inferred [O r1I]--H8 emission lines. No lensing correction has been applied. 


or displayed in the output maps. To produce the maps of the 
physical properties and study their spatial distribution, we 
display the 50th percentile of the inferred parameter, calculated 
from the posterior distribution that BAGPIPES provides. We can 
also present the uncertainties for each pixel extracted from the 
16th and 84th percentiles of the posterior distribution (see 
Carnall et al. 2018 for details). 


4. Results and Discussion 


In this section we present and discuss the results of our 
study. We provide both an integrated and a spatially resolved 
analysis of the physical properties that we infer for the five 
targets that compose this work. 


4.1. Spatially Resolved Physical Properties 


For all galaxies studied, we have produced maps of the 
physical parameters inferred with BAGPIPES. These include the 
SFR surface density (SFRD), the stellar mass surface density 
(SMD), the visual extinction Ay, the mass-weighted age, the 
UV slope 5, and the equivalent width (EW) of the [O m] and 
H8 emission lines. The last two are measured from the 
BAGPIPES posterior SEDs. We also display maps of the 
empirical colors F150W—F277W, which is a proxy for the UV 
slope, and F356W-F444W, which generally traces the strength 
of the inferred [O ri1]--H emission (except for the lowest- and 
highest-redshift galaxies, where no available couple of bands 
capture the lines and continuum accordingly). 

Figures 2—6 display the resulting maps for the five galaxies 
presented in this work, from the lowest redshift (z — 5.275) to 
the highest (z — 8.498). First, we see that all galaxies are 


resolved and display strong empirical color gradients, both in 
the blue bands and in the red bands. These gradients appear on 
larger scales than the FWHM of the F444W PSF (07145, 
equivalent to ~3.6 pixels) confirming that we can resolve 
trends and structures in our sample. In general, we find that, 
even at these early times, most of the galaxies display multiple 
star-forming clumps, traced by the F200W contours, as is 
found also by other recent works (e.g., Chen et al. 2023; 
Claeyssens et al. 2023; Treu et al. 2023). These are regions of 
very high inferred EWs of the [O M]+H8 emission (in the 
range ~300-4000 A rest frame), embedded within larger 
structures that are not undergoing a burst of star formation. In 
these regions with extreme line EWs, the inferred ages are 
extremely young («10 Myr), corresponding to a bursty clump 
of young stars that is resolved in targets ID10612 (Figure 4) 
and ID6355 (Figure 5) and marginally unresolved in ID5144 
(Figure 3) and ID4590 (Figure 6), given the scale of the F444W 
PSF FWHM. Around these high-EW bursty clumps, we also 
find underlying older stellar populations (~100 Myr), which 
would be missed in an integrated analysis, as we will discuss in 
more detail in the next subsection. 

Figure 2 displays the maps for galaxy ID8140, at redshift 
z= 5.275. We obtain smooth maps for all physical properties 
and colors. The SFRD and SMD appear entirely cospatial, and 
the UV slope traces perfectly the dust obscuration map. This 
galaxy has ~2 times older stellar populations compared to the 
average of the rest of the sample, in line with being the galaxy 
at lowest redshift. It shows strong color, Ay, and UV slope 
gradients, with two distinct clumps, one red and one blue. This 
could indicate that this galaxy is undergoing a merger, even 
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Figure 3. Maps of the empirical colors and the physical properties inferred with BAGPIPES on the galaxy ID5144 at z — 6.383. See Figure 2 for more details. 
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Figure 4. Maps of the empirical colors and the physical properties inferred with BAGPIPES on the galaxy ID10612 at z — 7.663. See Figure 2 for more details. 


though the EWs (the lowest of all targets) and ages show little 
variation across the object. 

Figure 5 shows the resulting maps for the galaxy ID6355 at 
z=7.665. It is the galaxy with the most extreme EWs 
(reaching ~4000 A). We see that the region with very high 
EW is very extended for this galaxy, leaving barely a shell 


where we find underlying older stellar populations, which are 
otherwise outshined (or not present) by the younger stars in 
these strong line emission regions. The clear gradients in the 
empirical colors give us confidence that this shell is real and 
not an artifact of the age—dust degeneracy in the SED fitting 
process, which we also discuss in further detail in Section 4.4. 
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Figure 5. Maps of the empirical colors and the physical properties inferred with BAGPIPES on the galaxy ID6355 at z — 7.665. See Figure 2 for more details. The three 
boxes A, B, and C indicate the pixels that are analyzed in more detail in the text, as well as in Figure 7, where the best-fit SEDs are shown for each individual pixel. 
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Figure 6. Maps of the empirical colors and the physical properties inferred with BAGPIPES on the galaxy ID4590 at z — 8.498. See Figure 2 for more details. 


On top of this, the shell is larger than the PSF scale. In Figure 7 the age and EW maps in Figure 5, since they appear to be very 


we present and analyze the fits for three individual pixels distinct regions within this galaxy. Albeit being toward the 
within this source, so that we can study further whether the edge of the galaxy, all three pixels fulfill our S/N threshold, so 
"shell" of older stars in this particular galaxy is real or an that we can produce robust fits (with reduced X? values within 
artifact. We select the pixels A, B, and C that are indicated in 0.10-0.53). These pixels are also far enough from each other so 
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Figure 7. Best-fit SEDs for the three pixels A, B, and C indicated in Figure 5, 
from the galaxy ID6355 at z — 7.665. The turquoise points and error bars 
correspond to the NIRCam photometry, the orange points to the best-fit model, 
and the orange curve and shaded region are the best-fit SED inferred with 
BAGPIPES and the corresponding 16th and 84th percentile uncertainty interval. 
The inferred physical parameters are indicated for each pixel, as well as the 
reduced x? of the fit. 


that the PSF is not blending the information they encode, and 
we can thus resolve their different stellar populations. We can 
clearly see that the SEDs look different, reflecting the gradients 
that we already see in Figure 5, both on the maps of the inferred 
physical parameters and on the empirical color maps. The 
greatest difference is observed in the strength of the inferred 
[O m]+H6 emission lines, since the SED for pixel C has 
extreme EW, reaching 4264 + 533 A rest frame. This yields a 
considerable difference in the inferred ages, with pixel A 
having a mass-weighted age of 159* 11? Myr and pixels B and C 
displaying very young stellar ages under 10 Myr. The corner 
plots for each fit can be found in Appendix C. 

Figure 6 shows the results for the galaxy ID4590, which is 
the highest redshift (z — 8.498) in this work. It is a very 
compact source, with only 31 pixels where S/N > 2 for both 
F150W and F200W. We see a star-forming clump, which also 
corresponds to the highest inferred stellar mass, and toward the 
dustiest zone in the Ay map. We see a marginally unresolved 
centrally located clump of young stellar population, which is 
not entirely cospatial with the star-forming burst but traces 
perfectly the higher EW of the inferred [O mM] and H8 emission 
lines. On top of this, the empirical color maps display a clear 
gradient, which follows the ones observed for the SFRD, SMD, 
Ay, and @ maps. The remaining targets exhibit similar trends 
for all physical properties, with the main characteristic being 
this region with extremely high EWs and therefore very young 
stellar populations. 


4.2. Integrated Analysis 


Besides providing an invaluable insight into the internal 
structure of galaxies, we want to test whether spatially resolved 
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Figure 8. Best-fit SEDs and models for the integrated (black curve and circles) 
and resolved (red curve and squares) modeling of the five galaxies studied in 
this work. The turquoise points and error bars correspond to the integrated 
NIRCam photometry. The red curve is inferred by summing the posterior 
distributions in all pixels and calculating the 50th percentile of the resulting 
one. The shaded regions correspond to the 16th and 84th percentiles of the 
summed posterior distribution. The inset cutouts correspond to the same RGB 
images as in Figure 1. The reduced X? values of each fit are indicated. 


observations yield other consequences, such as inferring 
different physical properties, compared to only integrated 
measurements. 

To perform this test, we sum the photometry in each 
Observed band for the pixels that fulfill our S/N criteria, so that 
we only consider the same pixels that we fit in the spatially 
resolved analysis shown in Figures 2-6. With the sum of the 
photometry in each filter, we then use BAGPIPES to find the 
best-fit SED and infer the integrated physical parameters. We 
use the exact same setup as for the spatially resolved run, 
described in Section 3.1. We present the integrated fits with the 
best-fit SEDs in Figure 8. The inferred integrated physical 
properties for each galaxy can be found in Table 3 in 
Appendix B. In Figure 8, we can see that both the resolved 
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Figure 9. Comparison between the stellar mass that we infer in the spatially 
resolved analysis and the integrated fit. No lensing correction is applied in any 
of the two estimates. The plotted values are shown in Table 2. The one-to-one 
line is indicated, as well as the 0.5 and 1 dex offset lines. 


and integrated models (red and black curves, respectively) fit 
adequately the photometry (turquoise points), with reduced x? 
values within the range 0.1—7.5. The surprising finding is that 
both best-fit SEDs are considerably different. For all galaxies, 
we find that the high EWs that we could spatially locate in the 
resolved analysis within a clump now completely dominate the 
overall fit. This results in inferring extremely young ages in the 
integrated light, and potentially too low stellar masses as a 
result. 

Figure 9 shows the comparison between the stellar mass 
estimates that we infer in the spatially resolved analysis and the 
integrated fit. Table 2 provides these values, as well as the 
mass-weighted ages. We obtain the spatially resolved mass 
estimate by summing the stellar mass inferred in each 
individual pixel. The resolved average age is inferred as the 
mean mass-weighted age of all pixels. In both cases, the 
resolved uncertainties reported in Table 2 are calculated with 
the 16th and 84th percentiles of the summed posterior 
distribution over all pixels. For the integrated run, the output 
from BAGPIPES is directly reported, and the uncertainties are 
calculated with the 16th and 84th percentiles from the posterior 
distribution. We find that the integrated run estimates system- 
atically lower stellar masses than the spatially resolved one, 
from ~0.5 up to | dex lower, seemingly without any trend with 
the redshift. This, as explained above, is a consequence of 
being forced to choose extremely young stellar populations to 
fit the integrated light, due to the strong emission lines 
dominating it completely. We see this in Table 2, since all 
galaxies have an integrated age of under 10 Myr, except the one 
at lowest redshift, with a best-fit age of 14.9*]j Myr. We 
demonstrate this by fixing the age in BAGPIPES to the average 
stellar age inferred in the spatially resolved analysis. By doing 
this, the estimate of the stellar mass in the integrated run 
increases, retrieving closer values to the spatially resolved 
masses. Moreover, we see that the average mass-weighted ages 
in the resolved analysis are all >10 Myr, and significantly older 
than the integrated ages, even when considering the large 
uncertainties associated with the age estimate. 
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Table 2 
Values for the Stellar Mass and Mass-weighted Stellar Age That We Infer with 
BAGPIPES, in Both the Integrated Run and the Spatially Resolved Analysis 
Presented in This Work 


3 ID log(uM. / Mc) Age (Myr) 
Integrated Resolved Integrated Resolved 
5.275 8140 9.147005 9.67*010 14.9*11 2511172 
6.383 5144 7.98071 8.907022 2.5*M 1477138 
7.663 10612 8.077 6:43 876031 1.6701 7148) 
7.665 6355 8.84005 9.19*023 1.3438 25*5 


8.498 4590 8.047041 9.09*022 1.9704 1114588 


Note. We plot the mass values in Figure 9. 


Figure 10 shows how the SFH affects the inferred stellar 
mass. We plot the sum of the SFH inferred for the spatial 
pixels, as well as the SFH estimated in the unresolved analysis, 
for the galaxy ID10612 at z=7.663. The integrated SFH 
consists of a single burst with very young age (~2 Myr), 
whereas the spatially resolved SFH is a distribution that covers 
a wider age range, reaching up to ~300 Myr. For this galaxy, 
this would mean a formation redshift of z~ 12. We see that 
whereas the integrated analysis forms all stellar mass within 
less than 10 Myr, this only corresponds to ~6% of the spatially 
resolved stellar mass, directly proving where the mass 
discrepancy is coming from. We obtain the same results in 
the SFH comparison for all galaxies studied in this work. 

This could considerably change our current picture of mass 
assembly in the early universe, particularly while our samples 
are limited to the brightest candidates in small numbers. These 
systematics would affect from the stellar mass functions that we 
have derived so far at high redshifts to our cosmological 
models of galaxy formation and mass buildup, since all our 
observations and mass estimates at high redshift until now have 
been based on integrated measurements. Overall, the nature of 
these galaxies can completely change by having a resolved 
picture. 

On top of this, we adopt a parametric constant SFH model. 
Some works have shown that using instead a nonparametric 
SFH model can lead to inferring larger stellar masses by up to 1 
dex (all in an integrated approach), in particular for galaxies 
like the ones studied here, where a burst with very young stellar 
ages dominates the integrated light (e.g., Leja et al. 2019; 
Lower et al. 2020; Tacchella et al. 2022; Whitler et al. 2023). 
Nonparametric SFHs therefore seem like the better option 
when considering integrated SEDs, since we have shown that 
there can be significant spatial variation in star formation (see 
Figure 10). Moreover, works like Bisigello et al. (2019) argue 
that we need to include the two MIRI bands in order to 
constrain better the stellar mass estimates in high-redshift 
galaxies, particularly in young galaxies with nebular emission 
lines, such as the ones we study here. 


4.3. Comparison with Other Works 


To put our results into context, we compare the physical 
parameters that we infer with other published works on these 
targets. As stated before, only integrated analyses are available 
in the literature. We expect our integrated measurements and 
estimates of some of the physical properties, such as the stellar 
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Figure 10. Comparison between the SFH that we infer in the spatially resolved 
analysis (red curve) and the integrated fit (black curve) for the galaxy ID10612 
at z — 7.663. The shaded regions correspond to the 16—84th percentile range in 
each case. 


mass and SFR, to be lower than any other work, since aperture 
photometry yields larger fluxes on all bands (~20% larger for 
the galaxy ID10612 at z= 7.663 if we use a 0"5 aperture), 
given that we only consider pixels that fulfill our S/N criteria. 
On top of this, previous works have used varying data 
reduction, zero-points, SED fitting codes, SFH models, 
attenuation curves, and magnification corrections, among other 
differences. Therefore, here we mostly focus on comparing the 
physical nature of the sources, as inferred by different works. 
The mass-weighted ages that we infer in our integrated 
analysis are all consistent within the uncertainties with those 
found for these same five targets by Carnall et al. (2023), from 
which we reproduce the SED fitting assumptions and 
parameters, except using a differently reduced photometry 
and a Calzetti et al. (2000) attenuation curve, as well as 
extending our nebular grid as explained in Section 3. 
Tacchella et al. (2022) study the stellar populations of three 
of our targets, the ones at highest redshifts z= 7.663, 
z—7.665, and z=8.498 (IDs 10612, 6355, and 4590, 
respectively). They use PROSPECTOR with a flexible nonpara- 
metric SFH prescription instead to model the SEDs. They find 
that the highest-redshift galaxy (ID4590) is undergoing a recent 
burst, inferring a young stellar age under 10 Myr, like we 
obtain in our integrated analysis. They cannot rule out older 
stellar ages in their analysis. In our spatially resolved maps (see 
the age map in Figure 6), we find that most pixels have a mass- 
weighted age above ~100 Myr, except the centrally located 
young burst, which dominates the integrated light. The 
particularly striking case is the z — 7.665 galaxy (ID6355). 
Tacchella et al. (2022) also infer an extremely young age of 
gs Myr, and they rule out the presence of older stellar 
populations, since the extreme [O rI1]--H lines dominate the 
emission. In this work we find a shell of older stars, confirmed 
by the empirical color gradients and the other tests discussed in 
Section 4.4 and Appendix A. Finally, for the target ID10612 at 
z = 1.663, they argue that its SFH, together with its morph- 
ology, could indicate that this source is undergoing a merger. 
This is consistent with what we find, in terms of both the 
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empirical color gradients and the clumpiness of the blue 
F150W and F200W bands. Moreover, we find two distinct 
populations within the galaxy in the spatially resolved maps 
(see Figure 4). Consistent with our results here, they find that 
inferring older stellar ages (in their case by adding emission- 
line constraints with NIRSpec spectra) leads to larger stellar 
masses of up to 1 dex for the galaxy ID4590 at z — 8.495, 
which is exactly what we find in our integrated versus spatially 
resolved comparison. 

Our work confirms the issue discussed in Tacchella et al. 
(2022), Topping et al. (2022), and Whitler et al. (2023), where 
young stars outshine and dominate the emission when 
compared to older stars, the presence of which is difficult to 
rule out via integrated measurements. This leads to inferring 
lower stellar masses. In Tacchella et al. (2022) they conclude 
that the SFH prior is of vital importance, and they only infer 
stellar ages older than 10 Myr in their fits when using a 
nonparametric SFH model with a continuous prior and older 
populations present. This leads to an increase in the stellar 
masses of up to 0.6 dex. Whitler et al. (2023) also use various 
SFH models to explore the potential presence of old stellar 
populations in seemingly young galaxies. They find stellar 
masses larger by up to an order of magnitude with nonpara- 
metric SFH versus constant SFH models. The impact of the 
SFH model in the inferred stellar mass has been studied by 
other works (see, e.g., Leja et al. 2019; Suess et al. 2022), with 
similar results. 

Outshining and its effects on stellar mass estimates have 
been studied at lower redshifts (see, e.g., Maraston et al. 2010; 
Pforr et al. 2012). Our results agree with previous works such 
as Sorba & Sawicki (2018), who find a discrepancy in the 
inferred stellar masses of up to a factor of ~5, when having 
resolved SED fitting, although their study only reaches z — 2.5. 
Moreover, they propose that unresolved studies should apply 
corrections to their mass estimates. This resolves the mass 
missing problem, in which a tension is found between the 
observed stellar mass density of the universe and the SFRD 
(see, e.g., Leja et al. 2015). By correcting stellar mass functions 
with resolved estimates, they find that these agree better with 
the observed star formation densities collected by Madau & 
Dickinson (2014). Leja et al. (2020) also solve this discrepancy 
by using flexible nonparametric SFH priors, which produce 
older ages, thus inferring ~50% higher stellar mass density. 

Endsley et al. (2022) study a population of UV -faint galaxies 
at a similar redshift range z ~ 6.5-8, also finding that the SEDs 
are dominated by young stellar populations, exhibiting low 
masses. They find the majority of their objects to appear very 
blue (8 ~ —2), with some dusty galaxies (8 ~ —1). With our 
integrated analysis, we find three galaxies with 8 ~ —2 and two 
targets with a value closer to —1. In the spatially resolved 
maps, we find values —2 < 8 < —1, with some very dusty 
regions reaching values around 3~ —0.5. 

At this redshift range, the majority of targets with high 
EW([Ornr]--H8) have an inferred young stellar age when 
considering a constant SFH, just as we find in our integrated 
analysis (see Figure 9 in Endsley et al. 2022). On the other 
hand, there are works that find evolved stellar populations 
(2100 Myr) at even higher redshifts, such as Furtak et al. (2023) 
with z~ 10-16 candidates in the SMACS 0723 field and 
Leethochawalit et al. (2023) with 7 «z «9 photometrically 
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selected galaxies in the GLASS-JWST ERS program, which 
infer a median mass-weighted age of 140 Myr. 

From a resolved point of view, in a sample of z~ 6-8 
galaxies in the Extended Groth Strip (EGS) field, Chen et al. 
(2023) find multiple clumps dominated by young stellar 
populations, as well as significant variations in the EW of the 
[O rir]--H lines. They find EWs with extreme values such as 
the ones we find for most of our targets (of the order 
7300-3000 A), which also yield young ages in their fits (as 
also found by Vanzella et al. 2023 in the Sunrise arc at z 6), 
confirming once more what we are finding for these high- 
redshift targets. Moreover, Pérez-González et al. (2023) also 
find strong [O rII]--H8 emission in a spatially resolved analysis 
using HST and JWST data from the CEERS survey in the EGS. 
They also link these findings with very young starbursts with 
possibly an underlying older stellar population. 

In summary, our results are consistent with the works that 
have been published so far studying these same galaxies, or 
targets at a similar redshift range with JWST. By integrating 
our maps, we can produce similar results and draw equivalent 
conclusions to the integrated works performed so far in these 
sources. By producing a spatially resolved analysis, we can 
demonstrate the presence of underlying older stellar popula- 
tions that are otherwise outshined in the integrated analyses, 
inferring larger stellar masses and considerably affecting our 
picture of the nature of these high-redshift galaxies. 


4.4. Caveats 


As briefly mentioned before, one could argue that the 
"shells" of older stellar populations where the young stars with 
high EWs are embedded could be instead a result of the dust 
—age degeneracy present in SED fitting software. To test 
whether the gradient in age is real, we perform a test in which 
we fix the extinction to the value given by Tacchella et al. 
(2022). We find a similar gradient, where there is a shell of 
older stellar populations surrounding the bursty young star- 
forming region. The effects of dust and age are now blended 
into the age map, since we fix the Ay, but the gradient persists. 
On top of that, the individual fits that we obtain fixing the dust 
obscuration are very poor, compared to leaving Ay as a free 
parameter in the BAGPIPES fit. This gives us confidence that our 
fits with Ay as a free parameter sample better the galaxy 
properties and that the gradient observed is real. In Appendix A 
we discuss in more detail the age uncertainties, focusing on the 
target ID6355 at z — 7.665. 

Another caveat that our spatially resolved analysis could 
have is whether the process of PSF matching affects our 
inferred maps. One could argue that the mass-weighted age 
map could result as an artifact of the PSF-matching procedure, 
where flux is redistributed radially to match the resolution of 
the F444W band. We only observe this radial distribution on 
the age and EW maps. We observe nonradial flat gradients 
across the galaxies on all of the remaining physical properties, 
as well as the empirical colors. These gradients extend across 
spatial scales larger than the FWHM of the F444W PSF. We 
therefore conclude that this is not an effect of PSF matching but 
a true young stellar population clump centrally located in most 
galaxies. 
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Finally, besides the S/N threshold that we impose in the 
noisiest bands, one could still doubt whether there is enough 
S/N per pixel to be able to infer robust physical parameters. To 
test this, we apply a Voronoi tesselation binning method on the 
targets, in order to achieve bins with a constant minimum S/N 
across the image and filters. Imposing a minimum S/N of 5 or 
even up to 10 in all bands, we find the same gradients and 
trends that we observe in the maps of the various inferred 
physical parameters in all galaxies. This, combined with the 
fact that our fits display good reduced x? values, gives us 
confidence that the S/N in each native pixel is sufficient to 
provide trustworthy estimates. 


5. Summary and Conclusions 


We present the first spatially resolved analysis of spectro- 
scopically confirmed 5 < z « 9 galaxies in the SMACS 0723 
ERO field. We use images in six bands obtained with NIRCam 
on board JWST, spanning the wavelength range 0.8-5 um. 
With the SED fitting software BAGPIPES, we model the SEDs 
on a pixel-by-pixel basis, being able to infer the physical 
parameters on a 180—240 pc scale. Our main findings and 
conclusions are the following: 


1. All galaxies are resolved and display strong empirical 
color gradients. Even at these early times, these galaxies 
display multiple star-forming clumps. 

2. We find regions that exhibit high EW of the [O m]+H6 
emission (up to ~3000-4000 À). These extreme star- 
bursts are embedded within regions with less specific star 
formation, which points to very bursty star formation 
happening on small scales («1 kpc), not galaxy-wide. 

3. The strong line emission regions dominate the integrated 
light, biasing the fits toward very young inferred ages of 
the stellar population («10 Myr). Only a resolved 
analysis demonstrates the presence of older stellar 
populations, which can be seen in the spatial maps. 

4. Resolving the stellar populations on a pixel-by-pixel 
basis leads to inferring from 0.5 up to ~1 dex larger 
stellar masses, when compared to an integrated analysis. 
Our analysis extends previous findings on the problem of 
outshining and its effects on stellar mass estimates, which 
so far has only been studied at lower redshifts (up 
to z^ 3). 


Current and upcoming observations with JWST will allow us 
to characterize the early universe and first galaxies in a new and 
more complete way. The combination of having confirmed 
redshifts with NIRSpec and the unprecedented resolution and 
depth of NIRCam imaging will transform how we study 
galaxies, changing our current views on their internal structure 
and mass assembly, among others. The systematics in stellar 
mass estimates found in this work would have strong 
implications in the shape and evolution of the stellar mass 
function at high redshift, particularly while samples are limited 
to small numbers of the brightest candidates. Furthermore, the 
process of galaxy formation could be more extended and earlier 
than previously thought, as is implied by the presence of 
evolved older stellar populations being outshone by the 
youngest stars. Only with a spatially resolved analysis can 
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we begin to untangle the complexity of the internal structure of 
galaxies at this epoch. 
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Appendix A 
Age Uncertainty 


As discussed in Section 4.4, the degeneracy between age and 
dust could yield uncertain estimates. This could be particularly 
concerning for the source ID6355 at z — 7.665, and one could 
argue that the shell that we see on the age map is not real. This 
would mean that there is no underlying older stellar 
population being outshined by the young stellar population 
that dominates the emission for this galaxy. Figure 11 
displays the 16th, 50th, and 84th percentiles of the mass- 
weighted age, obtained with the posterior distribution that 
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Figure 11. Mass-weighted age 16th, 50th, and 84th percentiles of the posterior 
distribution inferred with BAGPIPES on the galaxy ID6355 at z = 7.665. 


BAGPIPES infers for each individual pixel. In the 50th 
percentile, which is the one we choose to display for every 
inferred physical property in Figures 2—6, we see a shell of 
older stellar ages, as discussed in Section 4.1. With the 16th 
percentile image, we see that even if we assume the maximum 
lower uncertainty, the shell of old stars would still display 
ages above 10 Myr. That is still older than what is inferred by 
other works on this target, as well as in our integrated 
analysis, where we infer an age of 1.3*5$ Myr. With the 84th 
percentile image, we see that these old stars could be up to 
hundreds of megayears old. Therefore, even within the 
uncertainty range, we can confidently say that there are older 
stellar components present in this galaxy, opposite to what is 
concluded by Tacchella et al. (2022) and Carnall et al. (2023). 
This is only visible with a careful spatially resolved analysis. 
On the other hand, the very extended region, where the EW is 
extremely high, can only be fit by young stellar templates. 
Considering the uncertainties, we still only obtain young 
stellar populations in that region. 


Appendix B 
Integrated Properties 


Table 3 shows the resulting physical parameters inferred 
with BAGPIPES on the integrated fit for each galaxy. As a 
reminder, we expect these values to be lower than the ones 
inferred by other works, since we only consider the pixels 
here that fulfill a certain S/N criterion, instead of performing 
aperture photometry, which would yield larger fluxes and 
different physical estimates. The integrated run is performed 
like this to be able to do a one-to-one comparison with the 
spatially resolved analysis, as discussed in Section 4.2. 
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Table 3 
Resulting Physical Properties Inferred with BAGPIPES for the Integrated Measurements of the Five Galaxies 
Integrated Properties ID8140 ID5144 ID10612 ID6355 ID4590 
z7—5275 z = 6.383 z = 7.663 z= 7.665 z = 8.498 
log(M,./Mo) 9.14 055 798'$1 8.077008 8.84706 8.04.0, 
SFR („Mo yr )) 16.013 1.0303 1.2*03 6.9*03 L1*07 
Av (mag) 1.16002 0.587015 0.451044 1.017002 0.48008 
Age (Myr) 14.9*17 2.5114 1.6*07 1.3438 1.9*03 
EW([O rri]--H5) (A) 729 + 44 1996 + 386 3048 + 329 3884 + 252 2565 + 254 
UV slope 8 —1.20 + 0.03 —1.88 + 0.11 —2.08 + 0.09 —1.40 + 0.05 —2.03 + 0.09 
Appendix C respectively) of the galaxy ID6355 at z= 7.665. The pixels are 
Individual Fits shown in the maps of Figure 5, and the best-fit SEDs and 


physical properties are displayed in Figure 7. As we would 
expect, the fits for B and C are better constrained than for A, 
since the latter is the most outer pixel in the galaxy, thus with 


Here we present the corner plots obtained with BAGPIPES on 
the fits for the individual pixels A, B, and C (Figures 12-14, 
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Figure 12. Corner plot of the BAGPIPES fit on pixel A of the galaxy ID6355 at Figure 13. Corner plot of the BAGPIPES fit on pixel B of the galaxy ID6355 at 
z = 7.665. z = 7.665. 
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Figure 14. Corner plot of the BAGPIPES fit on pixel C of the galaxy ID6355 at 
z = 7.665. 


the lowest S/N, albeit fulfilling our criterion. For this pixel, the 
distributions of age, log(U), and metallicity are broad, but we 
can still constrain the stellar mass even at low S/N. 
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